# Подключение к Google drive
from google.colab import drive
drive.mount('/content/drive')
import numpy as np
import pandas as pd
# Загрузим набор данных
df = pd.read_csv('/content/drive/My Drive/Colab Notebooks/freMPL-R.csv', low_memory=False)
df = df.loc[df.Dataset.isin([5, 6, 7, 8, 9])]
df.drop('Dataset', axis=1, inplace=True)
df.dropna(axis=1, how='all', inplace=True)
df.drop_duplicates(inplace=True)
df.reset_index(drop=True, inplace=True)
df = pd.read_parquet('data/df.parquet.gzip', engine = 'fastparquet')
В предыдущем уроке мы заметили отрицательную величину убытка для некоторых наблюдений. Заметим, что для всех таких полисов переменная "ClaimInd" принимает только значение 0. Поэтому заменим все соответствующие значения "ClaimAmount" нулями.
NegClaimAmount = df.loc[df.ClaimAmount < 0, ['ClaimAmount','ClaimInd']]
print('Unique values of ClaimInd:', NegClaimAmount.ClaimInd.unique())
NegClaimAmount.head()
df.loc[df.ClaimAmount < 0, 'ClaimAmount'] = 0
Для моделирования частоты убытков сгенерируем показатель как сумму индикатора того, что убыток произошел ("ClaimInd") и количества заявленных убытков по различным видам ущерба за 4 предшествующих года ("ClaimNbResp", "ClaimNbNonResp", "ClaimNbParking", "ClaimNbFireTheft", "ClaimNbWindscreen").
В случаях, если соответствующая величина убытка равняется нулю, сгенерированную частоту также обнулим.
df['ClaimsCount'] = df.ClaimInd + df.ClaimNbResp + df.ClaimNbNonResp + df.ClaimNbParking + df.ClaimNbFireTheft + df.ClaimNbWindscreen
df.loc[df.ClaimAmount == 0, 'ClaimsCount'] = 0
df.drop(["ClaimNbResp", "ClaimNbNonResp", "ClaimNbParking", "ClaimNbFireTheft", "ClaimNbWindscreen"], axis=1, inplace=True)
pd.DataFrame(df.ClaimsCount.value_counts()).rename({'ClaimsCount': 'Policies'}, axis=1)
import plotly.express as px
fig = px.scatter(df, x='ClaimsCount', y='ClaimAmount', title='Зависимость между частотой и величиной убытков')
fig.show()
Для моделирования среднего убытка можем рассчитать его как отношение величины убытков к их частоте.
df.loc[df.ClaimsCount > 0, 'AvgClaim'] = df['ClaimAmount']/df['ClaimsCount']
Класс для общих случаев
class InsDataFrame:
''' Load data method '''
def load_pd(self, pd_dataframe):
self._df = pd_dataframe
''' Columns match method '''
def columns_match(self, match_from_to):
self._df.rename(columns=match_from_to, inplace=True)
''' Person data methods '''
# Gender
_gender_dict = {'Male':0, 'Female':1}
def transform_gender(self):
self._df['Gender'] = self._df['Gender'].map(self._gender_dict)
# Age
@staticmethod
def _age(age, age_max):
if pd.isnull(age):
age = None
elif age < 18:
age = None
elif age > age_max:
age = age_max
return age
def transform_age(self, age_max=70):
self._df['driver_minage'] = self._df['driver_minage'].apply(self._age, args=(age_max,))
# Age M/F
@staticmethod
def _age_gender(age_gender):
_age = age_gender[0]
_gender = age_gender[1]
if _gender == 0: #Male
_driver_minage_m = _age
_driver_minage_f = 18
elif _gender == 1: #Female
_driver_minage_m = 18
_driver_minage_f = _age
else:
_driver_minage_m = 18
_driver_minage_f = 18
return [_driver_minage_m, _driver_minage_f]
def transform_age_gender(self):
self._df['driver_minage_m'],self._df['driver_minage_f'] = zip(*self._df[['driver_minage','Gender']].apply(self._age_gender, axis=1).to_frame()[0])
# Experience
@staticmethod
def _exp(exp, exp_max):
if pd.isnull(exp):
exp = None
elif exp < 0:
exp = None
elif exp > exp_max:
exp = exp_max
return exp
def transform_exp(self, exp_max=52):
self._df['driver_minexp'] = self._df['driver_minexp'].apply(self._exp, args=(exp_max,))
''' Other data methods '''
def polynomizer(self, column, n=2):
if column in list(self._df.columns):
for i in range(2,n+1):
self._df[column+'_'+str(i)] = self._df[column]**i
def get_dummies(self, columns):
self._df = pd.get_dummies(self._df, columns=columns)
''' General methods '''
def info(self):
return self._df.info()
def head(self, columns, n=5):
return self._df.head(n)
def len(self):
return len(self._df)
def get_pd(self, columns):
return self._df[columns]
Создаем класс-наследник, в котором переопределяем некоторые методы, специфические для конкретной ситуации, и создаем новые
class InsDataFrame_Fr(InsDataFrame):
# Experience (weeks to years)
@staticmethod
def _exp(exp, exp_max):
if pd.isnull(exp):
exp = None
elif exp < 0:
exp = None
else:
exp * 7 // 365
if exp > exp_max:
exp = exp_max
return exp
# Marital status
_MariStat_dict = {'Other':0, 'Alone':1}
def transform_MariStat(self):
self._df['MariStat'] = self._df['MariStat'].map(self._MariStat_dict)
# Social category
def transform_SocioCateg(self):
self._df['SocioCateg'] = self._df['SocioCateg'].str.slice(0,4)
Подгружаем данные
data = InsDataFrame_Fr()
data.load_pd(df)
data.info()
Преобразовываем параметры
# Переименовываем
data.columns_match({'DrivAge':'driver_minage','LicAge':'driver_minexp'})
# Преобразовываем
data.transform_age()
data.transform_exp()
data.transform_gender()
data.transform_MariStat()
data.transform_SocioCateg()
# Пересечение пола и возраста, их квадраты
data.transform_age_gender()
data.polynomizer('driver_minage_m')
data.polynomizer('driver_minage_f')
Для переменных, содержащих более 2 значений:
NB: В H2O не рекомендуется использовать one-hot encoding. Тем не менее, используем здесь фиктивные переменные, чтобы в дальнейшем сохранить возможность сравнения результатов построенных моделей.
# Onehot encoding
data.get_dummies(['VehUsage','SocioCateg'])
data.info()
col_features = [
'driver_minexp',
'Gender',
'MariStat',
'HasKmLimit',
'BonusMalus',
'OutUseNb',
'RiskArea',
'driver_minage_m',
'driver_minage_f',
'driver_minage_m_2',
'driver_minage_f_2',
'VehUsage_Private',
'VehUsage_Private+trip to office',
'VehUsage_Professional',
'VehUsage_Professional run',
'SocioCateg_CSP1',
'SocioCateg_CSP2',
'SocioCateg_CSP3',
'SocioCateg_CSP4',
'SocioCateg_CSP5',
'SocioCateg_CSP6',
'SocioCateg_CSP7'
]
col_target = ['ClaimAmount', 'ClaimsCount', 'AvgClaim']
df_freq = data.get_pd(col_features+col_target)
df_ac = df_freq[df_freq['ClaimsCount'] > 0].reset_index().copy()
from sklearn.model_selection import train_test_split
# Разбиение датасета для частоты на train/val/test
x_train_c, x_test_c, y_train_c, y_test_c = train_test_split(df_freq[col_features], df_freq.ClaimsCount, test_size=0.3, random_state=1)
x_valid_c, x_test_c, y_valid_c, y_test_c = train_test_split(x_test_c, y_test_c, test_size=0.5, random_state=1)
# Разбиение датасета для среднего убытка на train/val/test
x_train_ac, x_test_ac, y_train_ac, y_test_ac = train_test_split(df_ac[col_features], df_ac.AvgClaim, test_size=0.3, random_state=1)
x_valid_ac, x_test_ac, y_valid_ac, y_test_ac = train_test_split(x_test_ac, y_test_ac, test_size=0.5, random_state=1)
Пусть $y$ – целевая переменная, $X$ – матрица объясняющих переменных, $\beta$ – вектор параметров модели.
Матрица $X$ составлена из всех векторов наблюдений $x_i$, каждый из которых представляет собой объясняющую переменную.
Основные компоненты обобщенной линейной модели:
Регрессия Пуассона обычно используется в случаях, когда зависимая переменная представляет собой счетные значения и ошибки предполагаются распределенными в соответствии с распределением Пуассона. Зависимая переменная должна быть неотрицательной.
Функции связи между таргетом и объясняющими переменными предполагается логарифмической: $$g(\eta) = \ln(\eta) \Rightarrow \hat {y} = e^{x{^T}\beta + {\beta_{0}}}.$$
Модель оценивается методом максимального правдоподобия, функция логарифма правдоподобия с учетом штрафа регуляризации эластичной сети имеет вид: $$\max_{\beta,\beta_0} \frac{1}{N} \sum_{i=1}^{N} \Big( y_i(x_{i}^{T}\beta + \beta_0) - e^{x{^T_i}\beta + {\beta_0}} \Big)- \lambda \Big( \alpha||\beta||_1 + \dfrac {1} {2}(1 - \alpha)||\beta||^2_2 \Big),$$
где
Тогда соответствующая метрика Deviance имеет вид: $$D = -2 \sum_{i=1}^{N} \big( y_i \text{ln}(y_i / \hat {y}_i) - (y_i - \hat {y}_i) \big).$$
NB: Ниже $\lambda$ не имеет отношения к вышеупомянутому одноименному параметру регуляризации.
Напомним, что функция вероятности для распределения Пуассона имеет вид: $$p(k;\lambda) = \frac{\lambda^k e^{-\lambda}}{k!}, \hspace{10pt} \lambda\in\mathbb{R}^{+}.$$ Также, для распределения Пуассона справедливо, что: $$\mathbb{E}\left[k\right] = Var(k) = \lambda.$$ Тогда, для оценивания коэффициентов нашей модели необходимо максимизировать правдоподобие (совместную условную вероятность при имеющихся данных), что данные имеют распределение Пуассона: $$p(y_1,\dots,y_n|x_1,\dots,x_n;\beta_0,\beta) = \prod_{i=1}^{N}\frac{e^{y_i(x_i{^T}\beta + {\beta_{0}})} e^{-e^{x_i{^T}\beta + {\beta_{0}}}}}{y_i!} = L(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n).$$ Для упрощения задачи оптимизации перейдем к логарифму правдоподобия: $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = \sum_{i=1}^{N}\left(y_i(x_i{^T}\beta + {\beta_{0}}) -e^{x_i{^T}\beta + {\beta_{0}}}-\ln(y_i!)\right).$$ Поскольку величина $\ln(y_i!)$ не зависит от выбора параметров, можно упростить задачу: $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = \sum_{i=1}^{N}\left(y_i(x_i{^T}\beta + {\beta_{0}}) -e^{x_i{^T}\beta + {\beta_{0}}}\right).$$ Далее решается задача оптимизации для определения параметров модели (пакеты реализуют оптимизацию численными методами): $$\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta_0} = 0,\\\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta} = 0.$$ Обычно минимизируется отрицательное правдоподобие, которое является выпуклой функцией.
Метрика Deviance представляет собой отношение правдоподобия между двумя моделями: рассматриваемой моделью и "идеальной" моделью, в которая бы идеально предсказывала бы зависимую переменную. $$Deviance = 2(\ell_{ideal} - \ell_{model})$$
В качестве такой "идеальной модели" может использоваться сама зависимая переменная. Тогда, логарифм правдоподобия "идеальной модели" для GLM с распределением Пуассона имеет вид: $$\ell_{ideal} = \sum_{i=1}^{N}\left(y_i \ln(y_i) -y_i-\ln(y_i!)\right).$$
Из приведенного выше вывода правдоподобия для рассматриваемой модели, мы можем записать, обозначив $\hat{y}_i = e^{x{^T}\beta + {\beta_{0}}}$: $$\ell_{model} = \sum_{i=1}^{N}\left(y_i \ln(\hat{y}_i) -\hat{y}_i-\ln(y_i!)\right).$$
Тогда получаем, $$Deviance = 2\sum_{i=1}^{N}\left(y_i \ln(y_i) -y_i - y_i \ln(\hat{y}_i) +\hat{y}_i\right) = -2\sum_{i=1}^{N}\left(y_i \ln(y_i/\hat{y}_i) - (y_i -\hat{y}_i)\right).$$
GLM с гамма-распределением используется для моделирования положительной непрерывной зависимой переменной, когда ее условная дисперсия увеличивается вместе со средним значением, но коэффициент вариации зависимой переменной предполагается постоянным.
Обычно GLM с гамма-распределением используются с логарифмической или обратной функциями связи: $$g(\eta) = \ln(\eta);\hspace{20pt}g(\eta) = \frac{1}{\eta}.$$
Модель оценивается методом максимального правдоподобия, функция логарифма правдоподобия (для обратной функции связи) с учетом штрафа регуляризации эластичной сети имеет вид: $$\max_{\beta,\beta_0} - \frac{1}{N} \sum_{i=1}^{N} \frac{y_i}{x{^T_i}\beta + \beta_0} + \text{ln} \big( x{^T_i}\beta + \beta_0 \big ) - \lambda \Big( \alpha||\beta||_1 + \dfrac {1} {2}(1 - \alpha)||\beta||^2_2 \Big),$$
где
Соответствующая метрика Deviance имеет вид: $$D = 2 \sum_{i=1}^{N} - \text{ln} \bigg (\dfrac {y_i} {\hat {y}_i} \bigg) + \dfrac {(y_i - \hat{y}_i)} {\hat {y}_i}.$$
NB: Ниже $\alpha$ и $\lambda$ не имеют отношения к вышеупомянутым одноименным параметрам регуляризации.
Напомним, что функция плотности вероятности для Гамма-распределения имеет вид: $$f(x;\alpha,\lambda) = \begin{cases} \frac{\alpha^\lambda}{\Gamma(\lambda)}x^{\lambda-1}e^{-\alpha x},&x \ge 0\\ 0,& x <0 \end{cases},\hspace{20pt} \Gamma(x) = \int_0^{\infty}x^{\lambda-1}e^{-x}dx. $$ Также, для Гамма-распределения справедливо, что: $$ \mathbb{E}[x] = \frac{\lambda}{\alpha},\hspace{15pt} Var(x) = \frac{\lambda}{\alpha^2}.$$
Для оценивания GLM удобно параметризовать данное распределение иначе: $$ \mathbb{E}[x] = \frac{\lambda}{\alpha} = \mu,\hspace{15pt} Var(x) = \frac{\lambda}{\alpha^2} = \frac{\mu^2}{\lambda}.$$ Тогда, $$f(x;\mu,\lambda) = \frac{1}{\Gamma(\lambda)}\alpha^\lambda x^{\lambda-1}e^{-\alpha x} = \frac{1}{\Gamma(\lambda)}\left(\frac{\lambda}{\mu}\right)^\lambda x^{\lambda-1}e^{-\left(\frac{\lambda}{\mu}\right)x} = \frac{1}{x\cdot\Gamma(\lambda)}\left(\frac{\lambda x}{\mu}\right)^\lambda e^{-\frac{\lambda x}{\mu}},\hspace{15pt} x\ge 0.$$
Также напомним, что для GLM с функцией связи $\ln(x)$ предсказание (оценка математического ожидания) имеет вид $\hat {y} = e^{x{^T}\beta + {\beta_{0}}}$.
Тогда, для оценивания коэффициентов нашей модели необходимо максимизировать правдоподобие (совместную условную вероятность при имеющихся данных), что данные имеют Гамма-распределение, в предположении известного параметра $\lambda$: $$p(y_1,\dots,y_n|x_1,\dots,x_n;\beta_0,\beta) = \prod_{i=1}^{N}\frac{1}{y_i\cdot\Gamma(\lambda)}\left(\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right)^\lambda e^{-\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}} = L(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n).$$ Для упрощения задачи оптимизации перейдем к логарифму правдоподобия: \begin{align} \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) &= \sum_{i=1}^{N}\left(-\ln(\Gamma(\lambda))-\ln(y_i)+ \lambda\ln\left(\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right) -\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right)\\ &= \sum_{i=1}^{N}\left(-\ln(\Gamma(\lambda))-\ln(y_i) + \lambda\ln\left(\lambda\right)+\lambda\ln\left(y_i\right)-\lambda\left({x_i{^T}\beta + {\beta_{0}}}\right) -\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right)\\ &= \sum_{i=1}^{N}\left((\lambda-1)\ln\left(y_i\right)-\ln(\Gamma(\lambda)) + \lambda\ln\left(\lambda\right)-\lambda\left({x_i{^T}\beta + {\beta_{0}}}\right) -\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right). \end{align} Тогда обозначим слагаемые, не зависящие от параметров модели ($\beta_0$, $\beta$), как некоторую константу $C$, получаем: $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = \sum_{i=1}^{N}\left(-\lambda\left(\left({x_i{^T}\beta + {\beta_{0}}}\right) +\frac{y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right) + C(y_i, \lambda)\right),$$ Таким образом, требуется максимизировать $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = -\sum_{i=1}^{N}\left(\left({x_i{^T}\beta + {\beta_{0}}}\right) + \frac{y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right),$$ Далее решается задача оптимизации для определения параметров модели (пакеты реализуют оптимизацию численными методами): $$\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta_0} = 0,\\\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta} = 0.$$ Обычно минимизируется отрицательное правдоподобие, которое является выпуклой функцией.
Метрика Deviance представляет собой отношение правдоподобия между двумя моделями: рассматриваемой моделью и "идеальной" моделью, в которая бы идеально предсказывала бы зависимую переменную. $$Deviance = 2(\ell_{ideal} - \ell_{model})$$
В качестве такой "идеальной модели" может использоваться сама зависимая переменная. Тогда, логарифм правдоподобия "идеальной модели" для GLM с распределением Пуассона имеет вид: $$\ell_{ideal} = -\sum_{i=1}^{N}\left(1 + \ln(y_i)\right).$$
Из приведенного выше вывода правдоподобия для рассматриваемой модели, мы можем записать, обозначив $\hat{y}_i = e^{x{^T}\beta + {\beta_{0}}}$: $$\ell_{model} = -\sum_{i=1}^{N}\left(\frac{y_i}{\hat{y}_i} + \ln(\hat{y}_i)\right).$$
Тогда получаем, $$Deviance = 2\sum_{i=1}^{N}\left(- 1 - \ln(y_i) + \frac{y_i}{\hat{y}_i} + \ln(\hat{y}_i)\right) = 2\sum_{i=1}^{N}\left(-\ln\left(\frac{y_i}{\hat{y}_i}\right)+\frac{y_i-\hat{y}_i}{\hat{y}_i} \right).$$
!apt-get install default-jre
!java -version
pip install h2o
import h2o
from h2o.estimators.glm import H2OGeneralizedLinearEstimator
h2o.init()
# Преобразование в H2O-Frame
h2o_train_c = h2o.H2OFrame(pd.concat([x_train_c, y_train_c], axis=1))
h2o_valid_c = h2o.H2OFrame(pd.concat([x_valid_c, y_valid_c], axis=1))
h2o_test_c = h2o.H2OFrame(pd.concat([x_test_c, y_test_c], axis=1))
# Инициализируем и обучим GLM модель c кросс-валидацией
glm_poisson = H2OGeneralizedLinearEstimator(family = "poisson", link = "Log", nfolds=5)
glm_poisson.train(y="ClaimsCount", x = h2o_train_c.names[1:-1], training_frame = h2o_train_c, validation_frame = h2o_valid_c)
# Параметры модели: распределение, функция связи, гиперпараметры регуляризации, количество использованных объясняющих переменных
glm_poisson.summary()
# Метрики качества модели - по всем данным и на кросс-валидации
glm_poisson.cross_validation_metrics_summary().as_data_frame()
# Таблица коэффициентов модели (в зависимости от модели могут выводиться также стандартная ошибка, z-score и p-value)
glm_poisson._model_json['output']['coefficients_table'].as_data_frame()
# Таблица нормированных коэффициентов по всем данным и на кросс-валидации
pmodels = {}
pmodels['overall'] = glm_poisson.coef_norm()
for x in range(len(glm_poisson.cross_validation_models())):
pmodels[x] = glm_poisson.cross_validation_models()[x].coef_norm()
coef = pd.DataFrame.from_dict(pmodels).round(5)
coef['overall'] = abs(coef['overall'])
coef.sort_values('overall',ascending=False)
# Построение прогнозных значений для обучающей, валидационной и тестовой выборок
c_train_pred = glm_poisson.predict(h2o_train_c).as_data_frame()
c_valid_pred = glm_poisson.predict(h2o_valid_c).as_data_frame()
c_test_pred = glm_poisson.predict(h2o_test_c).as_data_frame()
# Сохранение обученной модели
model_glm_poisson = h2o.save_model(model=glm_poisson, path="data/", force=True)
model_glm_poisson
# Преобразование в H2O-Frame
h2o_train_ac = h2o.H2OFrame(pd.concat([x_train_ac, y_train_ac], axis=1))
h2o_valid_ac = h2o.H2OFrame(pd.concat([x_valid_ac, y_valid_ac], axis=1))
h2o_test_ac = h2o.H2OFrame(pd.concat([x_test_ac, y_test_ac], axis=1))
h2o_train_ac
# Инициализируем и обучим GLM модель c кросс-валидацией
glm_gamma = H2OGeneralizedLinearEstimator(family = "gamma", link = "Log", nfolds=5)
glm_gamma.train(y="AvgClaim", x = h2o_train_ac.names[1:-1], training_frame = h2o_train_ac, validation_frame = h2o_valid_ac)
# Параметры модели: распределение, функция связи, гиперпараметры регуляризации, количество использованных объясняющих переменных
glm_gamma.summary()
# Метрики качества модели - по всем данным и на кросс-валидации
glm_gamma.cross_validation_metrics_summary().as_data_frame()
# Таблица коэффициентов модели (в зависимости от модели могут выводиться также стандартная ошибка, z-score и p-value)
glm_gamma._model_json['output']['coefficients_table'].as_data_frame()
# Таблица нормированных коэффициентов по всем данным и на кросс-валидации
pmodels = {}
pmodels['overall'] = glm_gamma.coef_norm()
for x in range(len(glm_gamma.cross_validation_models())):
pmodels[x] = glm_gamma.cross_validation_models()[x].coef_norm()
coef_ac = pd.DataFrame.from_dict(pmodels).round(5)
coef_ac['overall'] = abs(coef['overall'])
coef_ac.sort_values('overall',ascending=False)
# Построение прогнозных значений для обучающей, валидационной и тестовой выборок
ac_train_pred = glm_gamma.predict(h2o_train_ac).as_data_frame()
ac_valid_pred = glm_gamma.predict(h2o_valid_ac).as_data_frame()
ac_test_pred = glm_gamma.predict(h2o_test_ac).as_data_frame()
# Сохранение обученной модели
model_glm_gamma = h2o.save_model(model=glm_gamma, path="data/", force=True)
model_glm_gamma
h2o_df = h2o.H2OFrame(df_freq[col_features])
df['CountPredicted'] = glm_poisson.predict(h2o_df).as_data_frame()
df['AvgClaimPredicted'] = glm_gamma.predict(h2o_df).as_data_frame()
df['BurningCost'] = df.CountPredicted * df.AvgClaimPredicted
df[['CountPredicted', 'AvgClaimPredicted', 'BurningCost']].head()
# Разбиение датасета на train/val/test
df_ind = data.get_pd(col_features+['ClaimInd'])
x_train_ind, x_test_ind, y_train_ind, y_test_ind = train_test_split(df_ind[col_features], df_ind.ClaimInd, test_size=0.3, random_state=1)
x_valid_ind, x_test_ind, y_valid_ind, y_test_ind = train_test_split(x_test_ind, y_test_ind, test_size=0.5, random_state=1)
# Преобразование в H2O-Frame
h2o_train_c = h2o.H2OFrame(pd.concat([x_train_ind, y_train_ind], axis=1))
h2o_valid_c = h2o.H2OFrame(pd.concat([x_valid_ind, y_valid_ind], axis=1))
h2o_test_c = h2o.H2OFrame(pd.concat([x_test_ind, y_test_ind], axis=1))
# Преобразуем целевую переменную ClaimInd в категориальную при помощи метода asfactor во всех наборах данных
h2o_train_c['ClaimInd'].asfactor()
h2o_valid_c['ClaimInd'].asfactor()
h2o_test_c['ClaimInd'].asfactor()
h2o_train_c['ClaimInd']
# Инициализируем и обучим GLM модель c кросс-валидацией
glm_binomial = H2OGeneralizedLinearEstimator(family = "binomial", link = "Logit", nfolds=5)
# Инициализируем и обучим GLM модель c кросс-валидацией
glm_binomial.train(y="ClaimInd", x = h2o_train_c.names[1:-1], training_frame = h2o_train_c, validation_frame = h2o_valid_c)
# Параметры модели: распределение, функция связи, гиперпараметры регуляризации, количество использованных объясняющих переменных
glm_binomial.summary()
# Метрики качества модели - по всем данным и на кросс-валидации
glm_binomial.cross_validation_metrics_summary().as_data_frame()
# Таблица коэффициентов модели (в зависимости от модели могут выводиться также стандартная ошибка, z-score и p-value)
glm_binomial._model_json['output']['coefficients_table'].as_data_frame()
# Таблица нормированных коэффициентов по всем данным и на кросс-валидации
pmodels = {}
pmodels['overall'] = glm_binomial.coef_norm()
for x in range(len(glm_binomial.cross_validation_models())):
pmodels[x] = glm_binomial.cross_validation_models()[x].coef_norm()
coef_ac = pd.DataFrame.from_dict(pmodels).round(5)
coef_ac['overall'] = abs(coef['overall'])
coef_ac.sort_values('overall',ascending=False)
# Построение прогнозных значений для обучающей, валидационной и тестовой выборок
с_train_pred = glm_binomial.predict(h2o_train_c).as_data_frame()
с_valid_pred = glm_binomial.predict(h2o_valid_c).as_data_frame()
с_test_pred = glm_binomial.predict(h2o_test_c).as_data_frame()
Сохраним полученную модель
model_glm_binomial = h2o.save_model(model=glm_binomial, path="data/", force=True)
с_train_pred.iloc[:, 0]
from sklearn.metrics import accuracy_score, f1_score, confusion_matrix
# Выведем импортированные выше метрики классификации для обучающей, валидационной и тестовой выборок
cm_train_binominal = confusion_matrix(y_train_ind ,с_train_pred.iloc[:, 0])
cm_valid_binominal = confusion_matrix(y_valid_ind ,с_valid_pred.iloc[:, 0])
cm_test_binominal = confusion_matrix(y_test_ind ,с_test_pred.iloc[:, 0])
cm_valid_binominal, cm_train_binominal, cm_test_binominal
Как видно из матриц выше, у нас большое количество FN значений по всем данным.
f1_train_binominal = f1_score(y_train_ind ,с_train_pred.iloc[:, 0], average='weighted')
f1_valid_binominal = f1_score(y_valid_ind ,с_valid_pred.iloc[:, 0], average='weighted')
f1_test_binominal = f1_score(y_test_ind ,с_test_pred.iloc[:, 0], average='weighted')
print(f'train - {f1_train_binominal}, valid - {f1_valid_binominal}, test - {f1_test_binominal}')
accuracy_score
ac_train_binominal = accuracy_score(y_train_ind ,с_train_pred.iloc[:, 0])
ac_valid_binominal = accuracy_score(y_valid_ind ,с_valid_pred.iloc[:, 0])
ac_test_binominal = accuracy_score(y_test_ind ,с_test_pred.iloc[:, 0])
print(f'train - {ac_train_binominal}, valid - {ac_valid_binominal}, test - {ac_test_binominal}')
Метрики показывают, что модель получилось явно слабой и близка просто к угадыванию 50% на 50%. Думаю, что ситуацию можно исправить сгенерировать новые фичи. Так же выборка является не сбалансированной, что тоже плохо сказалось на модели.